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We develop a practical theoretical formalism for studying the critical properties of a trapped Bose-Einstein 
condensate using the projected Gross-Pitaevskii equation. We show that this approach allows us investigate the 
behavior of the correlation length, condensate mode and its number fluctuations about the critical point. Mo- 
tivated by recent experiments [Science 315, 1556 (2007)] we calculate the critical exponent for the correlation 
length, observe clear finite-size effects, and develop characteristic length scales to characterize the finite-size 
influences. We extend the Binder cumulant to the trapped system and discuss an experimental method for 
measuring number fluctuations. 
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I. INTRODUCTION 

Phase transitions normally arise out of the competition be- 
tween thermal fluctuations and inter-particle interactions. An 
important counterexample is that of the ideal Bose gas, where 
the second order transition to the condensed phase occurs in 
the absence of interactions, being driven solely by quantum 
statistical effects. However, the inclusion of even very weak 
interactions into the description of the condensation transition 
has proven to be a great theoretical challenge. For example, 
agreement on the effect of vanishingly small s-wave interac- 
tions on the critical temperature was only recently reached, 
after approximately 50 years of debate (e.g. see ID). 

Perhaps one of the most beautiful results of statistical 
physics is the concept of universality - that the properties of 
a system in the critical region are independent of the micro- 
scopic details of the system, and depend only on a few general 
system features such its dimensionality and the symmetry of 
the order parameter that emerges at the transition. In this con- 
text, the critical behavior of a weakly interacting 3D Bose gas 
should be identical to that of ''He at the superfluid transition, 
or to other systems of the same universality class (usually re- 
ferred to as the 3D XY universality class). Accordingly, high 
precision experimental measurements of critical ''He (e.g. see 
llJl) are usually compared against theoretical calculations us- 
ing idealized XY or 0^ models |3, 4, 5), rather than a micro- 
scopic description of the system. 

The availability of experimental techniques for measuring 
correlations |6, 7 8 9 10, 11 1 is an important feature of the 
ultra-cold atom systems that has received extensive theoretical 
attention, particularly in relation to (zero temperature) quan- 
tum phase transitions (e.g. see O [T3l fT4l fT5l fT6l [TTl fill [191 
l20l ). Another area of interest is the effect of critical fluctu- 
ations on correlations in the system at the finite temperature 
transition, particularly as other quantities usually examined in 
condensed matter (e.g. susceptibilities and heat capacity) are 
not easy to measure in atomic gases. While the critical fluctu- 
ations of weakly interacting ID and 2D Bose gases are domi- 
nant over a wide temperature range (e.g. see 12111221 ). it was 
previously thought that the width of the critical region about 
condensation temperature, Tc, would be far too narrow to per- 
mit experimental investigation in the 3D system. However, 
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Figure 1: Schematic view of the critical regime of a trapped Bose 
gas. The local chemical potential varies across the system due to the 
harmonic trap potential. When the system is critical at the center, the 
criticality extends over the finite spatial range L with a normal cloud 
boundary condition. 



in extraordinary recent experiments Donner et al. |23 1 have 
made such measurements of a trapped 3D Bose gas, and were 
able to determine the critical exponent for the divergence of 
the correlation length to be = 0.67 ± 0.13. In its own right 
this result is an impressive demonstration that universality ap- 
plies to a mesoscopic system with of order 10^ atoms. Ad- 
ditionally, this direct measurement of two-point correlations 
is of interest because it has not been possible to do this in 
helium, although a value for the correlation length critical ex- 
ponent of helium = 0.67056 ± 0.0006) is inferred from 
the heat capacity exponent using Josephson scaling relation. 

Finite-size effects have a profound influence on the critical 
properties of a system, and have been extensively studied to 
understand the cross-over of helium critical behavior during 
dimensional reduction |24|. Such systems, confined to a finite 
region of size L, are well-described by finite-size scaling the- 
ory Il25ll26l . This theory shows that there is a universal scaling 
function that relates physical quantities of the finite to infinite 
systems depending on the quantity, the ratio of the correlation 
length to the system size, and the nature of the boundary con- 
ditions. In this context the scenario occurring in the harmonic 
trap is rather interesting (see Fig. [T]), and was first considered 
by Damle et al. in 1996 |27 1. The effect of the trapping poten- 
tial is to slowly vary the local value of the chemical potential. 
If the gas is critical at trap center, then moving out radially. 
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generated and their properties are sampled. In Sec. Ill we 
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Figure 2: Schematic view of the c-field and the incoherent regions 
for a Bose gas in a harmonic trap potential, and the approximations 
we employ in our treatment of the collective mode dynamics. 



present the main results of our research. We describe the de- 
tailed macroscopic parameters of our simulations which span 
the critical point and present our results for the first order cor- 
relation functions and correlation length. In that section we 
also consider finite-size effects, and condensate number fluc- 
tuations. As far as we are aware, these are some of the first 
theoretical results systematically examining these quantities 
in the critical regime. Additionally, in the hope of stimulating 
further experimental work we present an observable, which 
we show is related to the condensate number fluctuations, and 
could be used to measure condensate statistics in experiment. 



We conclude the paper in Sec. IV 



II. FORMALISM 



the system gradually becomes normal. Thus the finite-size 
boundary conditions are rather difficult to describe, as they 
require an understanding of the (non-universal) normal sys- 
tem. The experiment by Dormer et al. [23 1 did not observe 
finite-size effects: their two-point measurements were made 
over a region much smaller than the spatial extent of the criti- 
cal region, and yielded a value of the critical exponent in line 
with the uniform system. 

In our opinion the study of the finite-size effects in this sys- 
tem should be rich and would be worth additional theoreti- 
cal and experimental investigation. However, a major imped- 
iment to the theoretical development is the lack of a practical 
formalism for studying the harmonically trapped Bose gas in 
the critical regime. In this paper we report on progress towards 
such a formalism. We develop the projected Gross-Pitaevskii 
equation (POPE) for application to the critical regime of the 
Bose-Einstein condensation transition. The PGPE method is 
a c-field technique |28| applicable to the study of finite tem- 
perature degenerate Bose gases. It includes interactions be- 
tween low energy modes of the gas non-perturbatively and is 
applicable in the critical region, e.g. see [29,, 30 , 31, 32 1. In- 
deed, PGPE predictions for the shifts in critical temperature 
ll29l [33l are in good agreement with other theoretical predic- 
tions in the uniform case [ 1 [ and experimental measurements 
in the trapped system [34[. While this formalism has success- 
fully predicted equilibrium properties for a degenerate Bose 
cloud, there has been little quantitative work at the critical 
point. While quantum Monte Carlo methods are now avail- 
able for the trapped Bose gas [35[, we expect that the sim- 
pler PGPE method to be far more efficient, easy to apply, and 
hence widely used. Indeed, the PGPE method is essentially an 
inhomogeneous theory. In this context the applicability of 
the PGPE to the critical regime is hardly surprising. However, 
unlike helium, where such models are thought to be inappli- 
cable except in the critical regime, the PGPE theory also pro- 
vides an accurate description of the (non-critical) condensed 
and normal behavior for the dilute Bose gas. 

The organisation of the paper is as follows. In Sec. [ll|we 
review the projected Gross-Pitaevskii formalism, discussing 
how equilibrium states of the finite temperature Bose gas are 



We take our system to be described by the second quantized 
Hamiltonian 



where 



(r), (1) 



(2) 



is the single particle Hamiltonian, ^'(r) is the quantum Bose 
field operator, and C/q = ATrh^a/m is the interaction strength, 
with a the s-wave scattering length. The trap potential is given 

as 
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Voir) = -m {u;lx^ + + ^Iz'^) 



(3) 



A. PGPE formalism 

We briefly outline the projected Gross-Pitaevskii equation 
(PGPE) formalism, which is developed in detail in Ref. [36[. 
The Bose field operator is split into two parts according to 



*(r) = Vc(r)+Vi(r), 



(4) 



where ?/;c is the coherent region c-field and ijji is the inco- 
herent field operator (see [28]). These fields are defined as 
the low and high energy projections of the full quantum field 
operator, separated by the energy ecut, as shown in Fig. |2] In 
our theory this cutoff is implemented in terms of the harmonic 
oscillator eigenstates {(/j„(r)} of the time-independent single 
particle Hamiltonian i.e. e„(/3„(r) = H^ifnir), with e„ the 
respective eigenvalue. The fields are thus defined by 



nec 

^i(r) = ^anipn{v), 



(5) 
(6) 
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where the a„ are Bose annihilation operators, the c„ are com- 
plex amplitudes, and the sets of quantum numbers defining the 
regions are 



C 
I 



{n 
{n 



En < Ecut}, 
> Ecut}- 



(7) 
(8) 



1. Choice ofC region 

In general, the applicability of the PGPE approach to de- 
scribing the finite temperature gas relies on an appropriate 
choice for ecut, so that the modes at the cutoff have an av- 
erage occupation of order unity. This choice means that the 
all the modes in C are appreciably occupied, justifying the 
classical field replacement a,i c„- In contrast the I region 
contains many sparsely occupied modes that are particle-like 
and would be poorly described using a classical field approxi- 
mation. Because our interest here is in the critical regime, ad- 
ditional care is needed in choosing C. Typically strong fluctu- 
ations occur in the infra-red modes up to the energy scale Uqti 
where n is the density. Above this energy scale the modes are 
well-described by mean-field theory (e.g. see the discussion 
in ll37l [38l ). For the results we present here, we choose to 
have 

Ccut ^ ksT > Uon. (9) 



B. Equilibrium states 

In this subsection we review our procedure for calculating 
finite temperature equilibrium properties of a trapped Bose 
gas. The basic approach is to treat the C and I regions as in- 
dependent systems in thermal and diffusive equilibrium. We 
discuss the treatment of these regions separately below. Fur- 
ther details on this procedure are given in Sec. 3 of |28 1. 



1. PGPE treatment ofC region 



The equation of motion for 4'c is the PGPE 



ih 



dt 



i?oV'c + 7'c{C/o|^c|Vc} 



where the projection operator 



^c{i^(r)} 



nGC 



^„(r) / dv' ^l{v')F{v'), 



(10) 



(11) 



formalises our basis set restriction of t^c to the C region. The 
main approximation used to arrive at the PGPE is to neglect 
dynamical couplings to the incoherent region |39|. 



An important feature of Eq. (10 1 is that it is ergodic, so that 



the microstates -00 evolves through in time form a sample of 
the equilibrium microstates, and time-averaging can be used 
to obtain macroscopic equilibrium properties. Our basic pro- 
cedure for finding equilibrium states consists of evolving the 



PGPE with three adjustable parameters: (i) the cutoff energy, 
ecut, that defines the division between C and I, and hence the 
number of modes in the C region; (ii) the number of C region 
atoms, A'c; (iii) the total energy of the C region, Eq. The last 
two quantities, defined as 



Nr 



dr|0c(r)|' 



(12) 



(13) 



are important because they represent constants of motion of 



the PGPE ( 10 1, and thus control the equilibrium state of the 
system. 

2. Randomized initial states 

Given choices of these parameters, randomized initial states 
are constructed satisfying those constraints, and are evolved 
according to Eq. ( 10 1. We now briefly describe how we pro- 



duce these initial states. 

We choose to make use of the Thomas-Fermi approxima- 
tion to the condensate mode 



77tf(x) 



I flTF - Vb(x) 
Un 



^(mtf- Vo(x)), (14) 



where 9{x) is the unit step function and 



~2 



Oho / 



2/5 



(15) 



is the Thomas-Fermi chemical potential f^Ol, with uj — 
(io^LOyiDzYl'^ and Oho = yjhjmuj. We can generate a state 
of desired energy by superimposing the Thomas-Fermi state 
with a (high energy) randomized state, ?7ran(r) (normalized to 
Nq), according to 



77£;(r) = po'7TF(r) +Pl'7ran(r), 



(16) 



where the variables olq and oi\ are adjusted to obtain the de- 
sired energy. In practice, ?7ia„ is approximately orthogonal to 
ryxF and we can take a\ — \J\ — |aop- We note that since 
the randomized field only spans the C region, it is most con- 
veniently formed by the construction indicated in Eq. (jSj), but 
with the c„ taken to be random complex numbers, globally 
scaled so that I]„ec l^nP = ^c- 

Typically evolution times of order 20 trap periods are 
used for the system to relax towards equilibrium [41J, before 
properties of the equilibrium states are sampled using time- 
averaging, as we discuss below. 



3. Obtaining equilibrium properties for the C region 

To characterize the equilibrium state in the C region it is 
necessary to determine the average density, condensate frac- 
tion, temperature and chemical potential. Many of these quan- 
tities are also important for characterizing the I region (see 
Sec.|nB4|. 
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The average density is obtained as a time average over the 
field microstates, i.e., 



"c(r) 



1 



Ms 



(17) 



where {tj} is a set of A/.; times (after the system has been 
allowed to relax to equilibrium) at which the field is sam- 
pled. We typically use ^ 7000 samples over the 140 trap 
periods of our simulation to perform such averages. The no- 
tation nc(r) emphasizes this is the density contribution from 
atoms residing in the C region, with the obvious property that 
iVc - / drnc(r). 

To find the mean condensate number, (A^cond) |50|, in our 
equilibrium state, we use the Penrose Onsager definition ll42l . 
that (iVcond) is given by the largest eigenvalue of the one-body 
density matrix 



I.e., 



J dr'GW(r,r')^e 



(#(r)7^c(r')), (18) 



Qd 

(r') = (Afco„d)?/'co„d(r). (19) 



Like the density, we can evaluate the one-body density ma- 
trix as a time-average. Because the one-body density matrix 
characterizes the emergence of the order parameter in the sys- 
tem at the critical point, we will examine G'^\r,r') further 
in Sec. IIIB to characterize correlations in the system at the 
transition. 

Derivatives of entropy, such as the temperature and chem- 
ical potential, are difficult to evaluate for interacting micro- 
canonical systems. In 1997 Rugh proved that a microcanoni- 
cal average (and hence time averaging) of appropriate quanti- 
ties constructed from the Hamiltonian could be used to cal- 
culate entropy derivatives P3l . The Rugh method applies 
to classical mechanical systems, such as the PGPE which is 
the equation of motion for a system descrbed by the classi- 



cal Hamiltonian (energy functional) ( 12 1. The detailed imple- 
mentation of the Rugh formalism for the PGPE is rather tech- 
nical and we refer the reader to Refs. 1,29. ,44] for additional 
details of this procedure. 



4. Meanfield treatment of the I region 

The average properties of the incoherent region can be cal- 
culated, at the level of Hatree-Fock meanfield theory, from the 
one-particle Wigner distribution 

^i^"^' P) = 4 T^l ' (20) 

h-^ exp(/3[eHF(r,p) - ^^\) - 1 

where 



eHF(r,p) - ;^- + Fo(r) + 2C/o(nc(r) + ni(r)),(21) 
2m 



is the Hartree-Fock energy, and /i is the chemical potential. In 
this semiclassical description r and p are treated as continu- 
ous variables. However, care needs to be taken to ensure that 



Eq. (20 1 is only applied to the appropriate region of phase 



space spanned by the incoherent region, i.e. single-particle 
modes of energy exceeding ecut- In phase space this region 
is 



f^i = -1 r,P : ^- V"o(r) > ecut 

2m 



(22) 



This allows us to calculate density the incoherent region atoms 



ni(r) 



dpM^i(r,p), 



(23) 



the number of incoherent region atoms A'l = J dr rii(r), and 
thus the total number of atoms in the system N — Nq + Ni. 

In Ref 1 301 it is shown the the full one-body density matrix 
for the system is given by 



G(i) (r , r') = Gg' (r , r') + G^'^ (r , r') , (24) 



where G^-* (r, r') is given in Eq. (18 i and 



Gg^(r,r') f dpe-'P^"-"'^^'' 
Jn, 



m(^,P 



(25) 



However, since we take Ecut > Uqh, then G^^\r, r') contains 
only normal system correlations that decay on the length scale 
of the thermal de Broglie wavelength, AdB — h/ ^/2TTmkBT, 
so that the contribution of G[^\r,r') is negligible for |r — 
r'l > AdB. 



ni. RESULTS 
A. Sampling equilibrium states across the transition region 

We now discuss our procedure for generating equilibrium 
states spanning the condensation transition. We fix the vari- 
ables Nq and ecut to define our system and the generate equi- 
librium states with various energy values (Eq) finely spaced 
over a range where the thermalized condensate fraction is of 
order 1%. Varying Eq in this way causes T to vary (as is de- 
sired), but also causes the total number of atoms to vary (see 
Fig.[3]|. 

For each simulation we calculate the temperature and to- 
tal atom number using the methods described in the previ- 
ous section. The results for these quantities for the case of a 
rubidium-87 system with LLjx,y — 27rx 129 Hz, ojz = 2ttx 364 
Hz, ecut = 'i2hiOx and Nq, = 7573 are shown in Figs. [3|a) 
and (b). We have chosen to use these parameters rather than 
those of the ETH experiment |23| (which was carried out in 
a weaker trap), because the higher critical temperature of our 
parameters allows us to better satisfy the validity conditions 
of the PGPE theory (|9|. 

For each energy we perform 20 simulations (using different 
random initial states) and the spread in results seen in Figs. 
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Figure 3: Macroscopic parameters for a critical trapped Bose gas. 
(a) Total number of atoms and (b) temperature as the PGPE en- 
ergy Ec is changed, (c) Condensate fraction as a function of the 
reduced temperature r/Tci (A'^). Parameters: Rubidium-87 system 
with uj^^y = 2tv X 129 Hz, cj^ = 27r x 364 Hz, e^ut ~ I8.9hujx and 
TVc = 7573. 



|3ja) and (b) for each energy is indicative of the typical un- 
certainties in the thermal parameters. These results also show 
that as we change Eq the total number of atoms in the system 
changes quite appreciably. It is therefore convenient to work 
in terms of the reduced temperature, T' = T/Td, where 



Tel = T, 



cO 



UJ 



1.33 iV5 ) T,o, (26) 

flho 



with 



OMhivN^/^, (27) 

u; = {ujxUjyLOzy/^, Co = {ujx + uly + LUz)/^, and Oho = 
^JKJmio, see ||45l . The two terms in brackets in Eq. (26 1 cor- 
respond to the finite-size (oc N^^l'^') and meanfield interac- 
tion (c>c iV^/^) shifts of the critical temperature, respectively. 

In Fig. |3jc) we show the condensate fraction from these 
simulations as a function of T' . We note that T' — \ does 
not identify the transition precisely enough for understanding 
critical properties, as the above expression for Td excludes 
meanfield effects beyond first order and does not account for 
any critical fluctuation effects. 




Figure 4: First order correlation function gc{^x) for, from highest 
to lowest curves, T' = {0.957, 0.962, 0.969, 0.981, 0.990}. Param- 
eters for calculations as in Fig. [5] Thermal de Broglie wavelength 
distance scale indicated for reference. 



B. Spatial correlations and the correlation length 

In this section we consider the behavior of spatial correla- 
tions and the correlation length across the transition region. 
The spatial inhomogeneity of the trapped system requires ex- 
plicit consideration. Unlike the homogeneous system, where 
spatial fluctuations are only a function of the separation be- 
tween coordinates, in the trapped case both coordinates are 
separately important. So to study the development of order in 
the trapped system, we choose to examine correlations sym- 
metrically about the trap center to minimize inhomogeneous 
effects |51 1. We do this by defining the normalized correlation 
function 



5c(Aa;) 



(28) 



where x is the unit vector in the x direction, and ric(O) is the 
C region density at trap center As discussed in Sec. |II B 4 



the spatial correlations over distances exceeding Ads are de- 
scribed completely by the C region one-body density matrix. 



We evaluate g'^^^ (r, r') by time-averaging (see Sec. 



IIB3I, 



but due to our system's symmetry in the xy-plane, we can 
improve the quality of our results for gc(Aa;) by making use 
of radial averaging. Examples of g^i^^x) are shown in Fig. |4] 
In the region of the phase transition the first order correlation 
function is expected to take the form 



5c (Ax) oc -^e 



-Ax/i 



(29) 



for Ax > AdB> where ^ is the correlation length of the sys- 
tem. The variation in the correlation length as the temperature 
approaches the critical value is given by 



(30) 



for the uniform system, where v is the relevant critcal expo- 
nent. 
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To determine ^ we fit (29i to our numerical results for 



gc{Ax) on the spatial range 0.5/im < Ax < 2.2fi m (i.e. 



1.5A, 



dB 



< Ax < 



6.4AdB)- This range matches that used 
by Donner et al. Il23l . and ensures that we avoid having 
to deal with normal correlations at small separations, and 
inhomogeneous/finite-size effects at larger separations (also 
see in Sec. |inB The values of ^ we obtain are shown 
in Fig. [s] where we see ^ growing rapidly as T' approaches 
^ 0.96 from above. At temperatures below this the quality of 
the fits used to determine ^ is quite poor and there is appre- 
ciable scatter in the data points for ^. This poor fit arises from 
the development of appreciable condensate in the system (e.g. 
see the coldest results shown in Fig. |4]l. In the uniform sys- 
tem the condensate is spatially uniform and is easily neglected 
in correlation functions, however in the trapped system it ap- 
pears at the transition point in a spatially localized mode with 
a size of order the oscillator length, which is difficult to distin- 
guish from the non-condensate correlations. We interpret the 
data for T' < 0.963 as being below the transition point and 
our f values extracted using fits to ( [29] l in this regime as being 
unreliable. 

We then fit expression ( [30| to the correlation length results 
with T^, V, and an overall constant of proportionality as fit- 
ting parameters. The fit for our data is shown in Fig. |5] with 
a value of v ^ 0.8 ± 0.12 and = 0.963. We notice that 
while the fit is reasonable, there appears to be a certain degree 
of rounding off in the divergence of ^ near the critical point. 
While our data has appreciable scatter (mainly due to uncer- 
tainty in T), we expect that this rounding off is primarily due 
to finite-size effects. Additionally, our large uncertainty in the 
critical exponent arises because of the difficulty in locating the 
precise value of Tc for our system. In principle the divergence 
of ^ marks Tc, however the rounding off of this divergence 
and the scatter in our results adds uncertainties to the precise 
value of Tc that is difficult to quantify without a theory for the 
finite-size behavior of ^. 

Our fit value for ly differs from that for the 3D XY model 
(i.e. « 0.67) which is expected to be of the same univer- 
sality class, but is within the error bars of the Donner et al. 
experiment, which reported = 0.67 ±0.13. The critical tem- 
perature identified by our fit (i.e. w 0.963) is also shifted 
downward from the prediction of ( [29| . A similar downward 
shift in Tc was found by Davis et al. [33| in their analysis 
of the trapped Bose gas, arising because meanfield effects are 



typically underestimated by the analytic expression ( 26 1 



1. finite-size effects 

As our above results motivate, an important issue to con- 
sider in the trapped system is the role of finite-size effects 
II27I . At fixed temperature, the Ginzburg criterion for the 
dominance of critical fluctuations requires that the chemical 
potential (/i) differs from the critical value (/Xc) by no more 
than 



S^=\fl- < 



16Tr^ma^k%T^f) 
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Figure 5: The behavior of the correlation length across the con- 
densation transition. Solid line is a fit using l |30^ with u = 0.8 and 
= 0.963. Other parameters as in Fig. H] Shaded region around the 
solid line indicates fits within the error bar range (1/ = 0.8 ± 0.12). 
Rectangular shaded region indicates points excluded as being below 
T^. Dashed lines show values of quantities discussed in the text. 



e.g. see fA5\. In the trapped system the effective system chem- 
ical potential varies spatially according to /i(r) = /i — Vo(r) 
and, taking /i(r = 0) = /ic, we can map the Ginzburg condi- 
tion to a spatial length scale over which the system is critical. 
This length scale (diameter), along direction xj — {x,y,z}. 



Lj ~ 8\/27ra— ^ — —. 

nujj 



(32) 



which we shall refer to as the Ginzburg length. This sets 
the maximum correlation length that can occur in the sys- 
tem, thus defining the relevant parameter for assessing finite- 
size effects, and takes this (its largest) value when the center 
of the system is at the critical point. For the case of two- 
point correlations (measured along the x-direction for defi- 
niteness) another important length scale in the trapped sys- 
tem is Aa;i„ax^ the maximum point separation (about trap 
center) used to measure the correlation length (i.e. the fit 
of exp{—Ax/£,)/Ax is made over the range AdB < Ax < 

AxYYiax.}- 



System 


AdB(rc) 




Lx 


V 


Expt. 


0.5 /im 


2.2 ^m 


20 Aim 


0.67 ±0.13 


Th. 1 


0.34 /xm 


2.2 pm 


9 Atm 


0.8 ±0.12 


Th. 2 


0.42 /im 


2.2 /im 


6 Aim 


0.8±0.12t 



Table I: A summary of the parameters for critical property measure- 
ments. Expt. values refer to those of Donner ef a/. |23|. Th. 1 refer 
to the values for the main theoretical resu lts prese nted in this paper 



and Th. 2 to the results presented in Sec. Ill B 1 
are not fit, see text for additional discussion. 



These numbers 



The following conditions are required to accurately mea- 
sure critical properties and minimise finite-size effects 



(31) 



AdB < AXmax < i^i 



(33) 
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Figure 6: The correlation lengtii for a smaller system with Nc = 
3.03 X 10^ and ecut ~ 23hijJx, with other parameters as described 
in Sec. |III A| Solid line and surrounding shaded region are not fits to 
the data, but are precisely the same as those used in Fig. |5](see text). 
Rectangular shaded region indicates points excluded as being below 
Tc- Dashed lines show values of quantities discussed in the text. 



The first inequaUty ensures that there is a reasonable distance 
scale over which correlation measurements can be made to ac- 
curately determine the correlation length. The second inequal- 
ity ensures that finite-size effects are minimized. Obviously, 
finite-size effects cannot be completely avoided, since the cor- 
relation length can never diverge in the finite system, and Lj 
sets the maximum value we might expect for ^. The values 
of these various quantities for experiments and our results are 
shown in Table U 



To examine the influence of finite-size effects, we have per- 
formed calculations for a system with the same trap parame- 
ters considered for the main results presented in this paper, but 
with fewer atoms. In this case the critical physics occurs at a 
temperature of 200 nK and a lower ecut value of 23hLLjx 
is used. The relevant parameters for this system are summa- 
rized as "Th. 2" in Table |l] revealing that for this system all 
three length scales in ( [33] l are similar Our results for the cor- 
relation length behavior of this system are given in Fig. [6] and 
show a rather striking broadening of the critical behavior, as 
compared to the previous case displayed in Fig. [5] As a re- 
sult, this data is difficult to fit to the infinite system result (|30]l 
for the purposes of extracting the critical exponent. Instead of 
fitting, we simply place the same curves used in Fig. |5](i.e. 
same T^, and error bars) on the data and observe that it pro- 
vides an acceptable characterization of these results also. For 
both cases (Figs. [5land[6| we see that fits to the normal diver- 
gent expression (30 1 are good for < Lx/2, but significantly 



depart from this fit for larger values of Since the values of 
Lx differ by roughly a factor of two between these calcula- 
tions, this suggests that Lj. is indeed the correct length scale 
for assessing finite-size effects. 
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Figure 7: Binder cumulant behavior across the condensation transi- 
tion. Parameters same as in Fig. [3] 



Condensate number fluctuations and the generalized 
Binder cumulant 



An important issue to deal with in the trapped system is 
the identification of the critical point, as this will be needed 
for a better understanding of the critical region and the future 
development of higher precision calculations. For the uruform 
Bose gas the transition point is conveniently identified using a 
Binder cumulant, defined as 



(34) 



where A^o the population of the zero-momentum (condensate) 
mode. This Binder cumulant characterizes condensate num- 
ber fluctuations, and takes the universal value of C"'* = 
1.2430 at the transition (see ||29l ). 

Here we propose a generalization of the Binder cumulant to 
the trapped system of the form 



/~ig — (-^cond) 
" " (A^cond)- 



(35) 



with A^cond the condensate mode occupation. Our procedure 
to analyze the condensate number fluctuations is as follows. 
The condensate (lowest energy normal mode), ?/'cond(r), is de- 
termined according to the Penrose-Onsager method described 
in Sec. |IIB 1 [ using the time-averaged density matrix. We then 
use this mode to determine the instantaneous condensate am- 
plitude by evaluating the inner product 



Q^cond 



(36) 



on every microstate used to sample system properties. We 
identify A'cond = |Q!cond(ij)P as the condensate number in 
this microstate, and by sampling over long times, we can ob- 
tain histograms of the condensate fluctuations. 

In Fig. |7]we show our results for across the critical re- 
gion. These exhibit a rather dramatic reduction in the number 
fluctuations from values approaching 0^ — 2 (expected for 
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Figure 8: Histograms of condensate/ground mode number fluctua- 
tions across the critical point. 



the normal system) at temperatures slightly above the transi- 
tion point, to values of ~ 1 below the transition. At our 
transition temperature of « 0.963 (as determined from fit- 
ting the critical exponent) we find a mean value of « 1.03, 
a value well-below the expected for the uniform system in the 
thermodynamic limit. This suggests that if the Binder cumu- 
lant is a useful quantity for characterizing the condensation 
transition in the trapped case, the critical value of is much 
lower than the uniform system. 

Investigations of the number fluctuations of the condensate 
across the transition are of interest in their own right, and may 
be suitable to the techniques available in ultra-cold atom ex- 
periments. The subject of condensate number fluctuations has 
been extensively discussed in the dilute gas BEC literature 
(e.g. see the review of |46| and references therein), partic- 
ularly the unphysical large fluctuations for the ideal gas pre- 
dicted within the grand canonical ensemble. 

In Fig. |8]we show histograms of the the condensate (low- 
est mode) number distribution across the transition. We see 
the development of coherence in the system as the shape of 
the distribution changes from being maximum at A'cond — 
above the transition, to having a maximum at finite iVcond be- 
low the transition. We also find that as the temperature de- 
creases the condensate number fluctuations are suppressed, 
i.e. 



cond 



(iVcond))') 



0, 



(37) 



cond/ 



(as was implicit in the behavior of observed above) and 
that the distribution has negative skew. 

Given the phenomenal recent interest in measuring spatial 
correlations in ultra-cold atom systems |l6l |7] [8] |9] [TOl [TD. 
it would be of great interest to develop analogous techniques 
for observing these condensate number distribution in exper- 
iments. It is difficult to devise an experimental procedure 
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Figure 9: Comparison of (a-e) histograms of condensate/ground 
mode number distribution and (f-j) the the central momen- 
tum column density distribution across the critical point. Re- 
sults (a-e) and (f-j) correspond to the temperatures T' = 
{0.957, 0.962, 0.969, 0.981, 0.990} respectively. 



which could be used to measure acond (or A'cond) in a manner 
equivalent to ( [36] l, which requires complete phase and ampli- 
tude information about the field. So here we propose a quan- 
tity that can be directly measured in experiments and used 
to reveal the transition from incoherent to coherent number 
statistics of the condensate mode. In particular, we consider 
the central momentum column density 



dp 2 n(p) 



(38) 



as an observable, since it proportional to the peak density mea- 
sured in the usual absorption images taken of ultra-cold sys- 
tems WT\ . The motivation for choosing this quantity is that the 
long range coherence of the condensate is clearly revealed as 
a peak in momentum space, thus the central momentum value 
is correlated with the condensate occupation. The detailed re- 
lationship between nj,=o and A'cond is not unique, due to the 
contribution of the noncondensate to rtp=o- So measurements 
of rtp=o cannot be considered equivalent to the condensate 
yet, as we show below, there are clear qualitative similarities 
between the distributions of both quantities. 

In Fig. [9]we compare the distributions for rip^o and iVcond 
obtained from analysis of data sets from the same PGPE cal- 
culations. Qualitatively, the behavior of these distributions 
appears to be quite similar. The np=o distribution is clearly 
seen to be offset from zero at high temperatures as compared 
to the A'cond distribution (see Figs. [9|e) and {])). This off- 
set is related to the average momentum column density of the 
noncondensate component. 
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IV. CONCLUSIONS AND OUTLOOK 

In this work we have developed a theoretical tool for study- 
ing critical physics in the 3D trapped Bose gas, motivated by 
the recent theoretical work by the ETH group |23|. For nu- 
merical convenience we have studied a trapped Bose gas with 
different parameters to the experimental system, yet obtain a 
critical exponent that agrees to within the error bars of both re- 
sults. We have discussed finite-size effects and show that they 
can significantly alter the critical physics for systems with a 
small Ginzburg length. Finally, we considered fluctuations of 
the condensate mode occupation across the transition region, 
and have shown that measurements of the central momentum 
column density can be used to experimentally reveal the emer- 
gence of coherent statistics in the system. 

Our study here represents the first quantitative theoretical 
calculations for this system beyond meanfield level. However, 
there are numerous avenues for improvement of our approach 
that could be used to obtain higher precision results, which 
we believe will be needed to develop a better understanding of 
the trapped system critical physics, and stimulate more experi- 
ments in this fascinating new area for ultra-cold atomic gases. 
For instance, it would be interesting to apply the stochastic 
PGPE (SPGPE) formalism f48l to the critical regime, as both 



the chemical potential and temperature are control parameters 
in this approach. This may allow more precise measurement 
of critical exponents and would make direct comparisons with 
the fixed N experiment of Donner et al. |23| more easy to 
achieve I.52J . Another avenue of investigation would be to 
build on our formalism a more efficient method of sampling, 
e.g. using Monte Carlo algorithms (e.g. see |'37^). With ad- 
ditional improvements in precision we believe our formalism 
will be able to provide a detailed characterization of the finite- 
size cross over functions for the trapped Bose gas. Knowledge 
of these functions would be useful for several problems of cur- 
rent interest in the ultra-cold atomic physics, such as better 
understanding of the quasi-2D behavior and the emergence of 
phase defects in quenches across the critical regime [49 1 . 
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